clear all;
clc;

% ------------------------- File Information ------------------------------

numpat = 5;

epitag = 'neg_epi_bp_';
subtag = 'neg_sub_bp_';

epi_f  = [0.3,0.45,0.23,0.21,0.25];
sub1_f = [0.2,0.12,0.25,0.3,0.16];
sub2_f = [0.8,0.88,0.75,0.7,0.84];

noaxons = 200;
th_val  = 0;

% ------------------------- Importing Data --------------------------------

pol = -1;

sel_epi  = zeros(numpat,noaxons);
sel_sub1 = zeros(numpat,noaxons);
sel_sub2 = zeros(numpat,noaxons);

for i = 1:numpat
    
    ptag=['p',num2str(i),'_'];
    
    cd(['patient',num2str(i),'_thresholds']);

    epiloc  = ['_d',num2str(epi_f(i)),'_t',num2str(th_val)];
    epiloc  = strrep(epiloc,'.','p');
    epiloc  = strrep(epiloc,'-','n');
    sub1loc = ['_d',num2str(sub1_f(i)),'_t',num2str(th_val)];
    sub1loc = strrep(sub1loc,'.','p');
    sub1loc = strrep(sub1loc,'-','n');    
    sub2loc = ['_d',num2str(sub2_f(i)),'_t',num2str(th_val)];
    sub2loc = strrep(sub2loc,'.','p');
    sub2loc = strrep(sub2loc,'-','n');    
    
    epi_dc = [ptag,epitag,'dc_act',epiloc,'.txt'];
    epi_dr = [ptag,epitag,'dr_act',epiloc,'.txt'];
    epi_vth_dc = load(epi_dc);
    epi_vth_dr = load(epi_dr);
    epi_vth_dc = sort(pol*epi_vth_dc(:,1));
    epi_vth_dr = sort(pol*epi_vth_dr(:,1));
    
    sub1_dc = [ptag,subtag,'dc_act',sub1loc,'.txt'];
    sub1_dr = [ptag,subtag,'dr_act',sub1loc,'.txt'];
    sub1_vth_dc = load(sub1_dc);
    sub1_vth_dr = load(sub1_dr);
    sub1_vth_dc = sort(pol*sub1_vth_dc(:,1));
    sub1_vth_dr = sort(pol*sub1_vth_dr(:,1));
    
    sub2_dc = [ptag,subtag,'dc_act',sub2loc,'.txt'];
    sub2_dr = [ptag,subtag,'dr_act',sub2loc,'.txt'];
    sub2_vth_dc = load(sub2_dc);
    sub2_vth_dr = load(sub2_dr);
    sub2_vth_dc = sort(pol*sub2_vth_dc(:,1));
    sub2_vth_dr = sort(pol*sub2_vth_dr(:,1));

    cd ..;
    
    epi_dcvsdr  = pera_vs_perb(epi_vth_dc,epi_vth_dr);
    sub1_dcvsdr = pera_vs_perb(sub1_vth_dc,sub1_vth_dr);
    sub2_dcvsdr = pera_vs_perb(sub2_vth_dc,sub2_vth_dr);
    
    sel_epi(i,:)  = epi_dcvsdr(:,2)';
    sel_sub1(i,:) = sub1_dcvsdr(:,2)';
    sel_sub2(i,:) = sub2_dcvsdr(:,2)';  

end

sel_epi_mu = mean(sel_epi(5,:));
sel_epi_sd = std(sel_epi);

sel_sub1_mu = mean(sel_sub1(5,:));
sel_sub1_sd = std(sel_sub1);

sel_sub2_mu = mean(sel_sub2(5,:));
sel_sub2_sd = std(sel_sub2);

sel_epi_p=sel_epi(numpat,:);
sel_sub1_p=sel_sub1(numpat,:);
sel_sub2_p=sel_sub2(numpat,:);

% ------------------------- Plotting --------------------------------------

msub = [1,2:4:noaxons,noaxons];

f_dc = (100/noaxons)*(1:1:noaxons);

figure;
hold on;
plot(0:100,0:100,'k--','LineWidth',4);
plot(sel_epi_p(msub),f_dc(msub),'k-','LineWidth',4,'MarkerSize',10);
plot(sel_sub1_p(msub),f_dc(msub),'ko-','LineWidth',2,'MarkerSize',10,...
     'MarkerFaceColor','k');
plot(sel_sub2_p(msub),f_dc(msub),'ko-','LineWidth',2,'MarkerSize',10);
hold off;
% ttl = ['Selective Activation of Dorsal Columns (\theta = ',num2str(th_val),'^o)'];
% title(ttl,'FontSize',30,'FontWeight','b');
ylabel('Dorsal Columns Activated (%)','FontSize',36);
xlabel('Dorsal Roots Activated (%)','FontSize',36);
L1 = '1 mm above dura';
L2 = '1 mm above cord';
L3 = '1 mm below dura';
legend('No Selectivity',L1,L2,L3,'Location','NW');
set(gca,'FontSize',36);
set(gca,'XTick',0:25:100,'YTick',0:25:100);
axis tight;
axis square;
